clear
close all
clc

addpath(genpath('Estimation\lmjtoolbox'))

% Output folder for figures
Figures_path=strcat('.\Figures');

% Graphic settings
set(0,'defaultlinelinewidth',2.5);
set(0,'defaultaxesfontsize',14);  %12
set(0,'defaulttextfontsize',12);  %12
set(0,'defaultlinelinewidth',2.5);
rrs = get(0,'screensize');
rrf = get(0,'defaultfigurepos');
defpos = [00, 00, 1200, 500]; %2100 700
set(0,'defaultfigurepos',defpos);

% Settings
long_sample=1;  % 0 Sample ends in 2020:Q4 with vintage data(for the out-of-sample forecasts - Figure IX); 
                % 1 Estimation up to 2022:Q3 -> Figure IV, V, VI, VII, VIII
nirf=40;        % Time horizon impulse response functions (Figure IV)

%% Upload Quantitative model estimation output
if long_sample==0
    %% Model in Sample 1
    load ('.\Estimation\BFM_mod\ver01\FirstSampleOutput\end 2007.75 csminwel BFMPrior 05 09 23\workspace');
    KFStruS1=KFStru;
    [GG1,RR1,Cobs1,eu1,SDX1,ZZ1,initss1,ssvec1,flag1,ssnames1,namesstates1,namesshocks1,~,~,~,~] = ...
        modBFM_modver01(model.param,0,addsol);
    
    %% Model in Sample 2
    load('.\Estimation\BFM_mod\ver01\SecondSampleOutput\end 2020.75 csminwel BFMPrior 05 09 23\workspace.mat')
    KFStruS2=KFStru; 
    [GG2,RR2,Cobs2,eu2,SDX2,ZZ2,initss2,ssvec2,flag2,ssnames2,namesstates2,namesshocks2,~,~,~,~] = ...
        modBFM_modver01(model.param,0,addsol);
    
else %long sample
    %% Model in Sample 1
    load ('.\Estimation\BFM_mod\ver01\FirstSampleOutput\end 2007.75 csminwel BFMPrior 02 03 23\workspace');
    KFStruS1=KFStru;
    [GG1,RR1,Cobs1,eu1,SDX1,ZZ1,initss1,ssvec1,flag1,ssnames1,namesstates1,namesshocks1,~,~,~,~] = ...
        modBFM_modver01(model.param,0,addsol);
    
    %% Model in Sample 2
    load('.\Estimation\BFM_mod\ver01\SecondSampleOutput\end 2022.50 csminwel BFMPrior 02 03 23\workspace.mat');
    KFStruS2=KFStru;
    [GG2,RR2,Cobs2,eu2,SDX2,ZZ2,initss2,ssvec2,flag2,ssnames2,namesstates2,namesshocks2,~,~,~,~] = ...
        modBFM_modver01(model.param,0,addsol);
end



%% Figure IV
% Impulse Response Functions 
imp=zeros(size(RR2,2),1); imp2=zeros(size(RR2,2),1); imp3=zeros(size(RR2,2),1);
% Funded shock position 
imp(2)=1; 
%Unfunded shock position
imp2(3)=1; 
% Persistent cost-push shock position 
imp3(11)=1; 
%Initialization
ssirf=zeros(size(GG2,1),nirf); yyirf=zeros(size(ZZ2,1),nirf); ssirf2=zeros(size(GG2,1),nirf); yyirf2=zeros(size(ZZ2,1),nirf); 
ssirf3=zeros(size(GG2,1),nirf); yyirf3=zeros(size(ZZ2,1),nirf);
for i=0:nirf-1
%     Funded Transfers shock
    ssirf(:,i+1)=GG2^i*RR2*SDX2*imp; 
    yyirf(:,i+1)=ZZ2*ssirf(:,i+1);
%     UnFunded Transfers shock
    ssirf2(:,i+1)=GG2^i*RR2*SDX2*imp2;
    yyirf2(:,i+1)=ZZ2*ssirf2(:,i+1);
%     Persistent cost-push shock
    ssirf3(:,i+1)=GG2^i*RR2*SDX2*imp3;
    yyirf3(:,i+1)=ZZ2*ssirf3(:,i+1);
end
% Transpose the matrices containting the impulse responses to the shocks
ssirf=ssirf'; ssirf2=ssirf2';ssirf3=ssirf3';
yyirf=yyirf'; yyirf2=yyirf2'; yyirf3=yyirf3';

if long_sample==1
    fig=figure;
    
    hax=subplot(1,3,1);
    plot(1:nirf,4*yyirf(:,7),'--k')
    hold on
    plot(1:nirf,4*yyirf2(:,7),'b')
    hold on
    plot(1:nirf,4*yyirf3(:,7),'-.r')
    hold on
    plot(1:nirf,4*yyirf2(:,7)*0,'--m','LineWidth',.5)
    hold off
    axis tight
    title('Inflation')
    h1=legend('Funded Transfers','Unfunded Transfers','Persistent cost-push shock');
    set(h1,'FontSize',10,'Location','NorthEast')
    ax=subplot(1,3,2);
    
    plot(1:nirf,ssirf(:,20)*4,'--k')
    hold on
    plot(1:nirf,ssirf2(:,20)*4,'b')
    hold on
    plot(1:nirf,ssirf3(:,20)*4,'-.r')
    hold on
    plot(1:nirf,4*yyirf2(:,7)*0,'--m','LineWidth',.5)
    hold off
    axis tight
    title('Real Interest Rate')
    
    hax(3)=subplot(1,3,3);
    plot(1:nirf,yyirf(:,10)/4,'--k')
    hold on
    plot(1:nirf,yyirf2(:,10)/4,'b')
    hold on
    plot(1:nirf,yyirf3(:,10)/4,'-.r')
    hold on
    plot(1:nirf,4*yyirf2(:,7)*0,'--m','LineWidth',.5)
    hold off
    axis tight
    title('Debt-GDP Ratio')
    savefigure_pdf(strcat(Figures_path,'\Figure_IV'));
end
%% Upload the Data Set

if long_sample==0
    Dataf=xlsread('.\Estimation\BFM_mod\data\usdatav1');  % Real (older vintage) data to replicate Figure IX.
    Data=Dataf(1:end-1,:); %delete the last row (2021:Q1 observation, which will be useful for later)
else 
    Data=xlsread('.\Estimation\BFM_mod\data\usdata');  
end

TT=Data(4:end,1)'; YY=Data(4:end,2:end);
nobs=length(TT);
msel=[ones(192,1);ones(nobs-192,1)*2]; %Dummies capturing the change in the estimated parameters from the first and the second sample

%% Kalman Smoother
DD=Cobs1([144,70:71,7,3,73,14,74,72,145,3*ones(1,10),14]); 
ZZ_t=zeros(size(ZZ2,1),1); HH=zeros(size(ZZ2,1),size(ZZ2,1));

GG(:,:,1)=GG1;GG(:,:,2)=GG2;
RR(:,:,1)=RR1;RR(:,:,2)=RR2;
QQ(:,:,1)=SDX1.^2;QQ(:,:,2)=SDX2.^2;

[ SSt,AAt,PPt] = KalmanSmoother( YY,GG,RR,QQ,DD,ZZ2,HH,ZZ_t,msel );

SSt=SSt';AAt=AAt';

%% Extract the Smoothed Innovations

% Append Shocks
GGi(:,:,1)=[[GG1,zeros(size(GG1,1),size(RR1,2))];[zeros(size(RR1,2),size(GG1,2)+size(RR1,2))]];
GGi(:,:,2)=[[GG2,zeros(size(GG2,1),size(RR2,2))];[zeros(size(RR2,2),size(GG2,2)+size(RR2,2))]];
RRi(:,:,1)=[RR1;eye(size(RR1,2))];
RRi(:,:,2)=[RR2;eye(size(RR2,2))];
YYi=SSt;
ZZi=[eye(length(GG2)),zeros(size(GG2,1),size(RR2,2))];DDi=zeros(length(GG2),1); ZZ_ti=DDi;
Ati=zeros(length(GGi),1); Pti=eye(length(GGi));
HHi=eye(size(ZZi,1))*.0001;
[ SSti,AAti,PPti] = KalmanSmoother( YYi,GGi,RRi,QQ,DDi,ZZi,HHi,ZZ_ti,msel,Ati,Pti );
SSti=SSti';

%% Simulating the Model One Shock at the Time (various counterfactual simulations)

sst=SSt(1,:)'; ssAt=sst; ssBt=sst; ssCt=sst; ssDt=sst; ssICt=sst;
for t=2:nobs
    sst(:,t) =squeeze(GG(:,:,msel(t)))*sst(:,t-1)+squeeze(RR(:,:,msel(t)))*SSti(t,size(SSt,2)+1:end)'; %Simulation with all the estimates shocks
    ssAt(:,t)=squeeze(GG(:,:,msel(t)))*ssAt(:,t-1)+squeeze(RR(:,3,msel(t)))*SSti(t,size(SSt,2)+3)'; %Simulation with only the unfunded transfers shocks
    ssBt(:,t)=squeeze(GG(:,:,msel(t)))*ssBt(:,t-1)+squeeze(RR(:,[4:5,7:11],msel(t)))*SSti(t,size(SSt,2)+[4:5,7:11])'; %Simulation using only non-policy shocks
    ssICt(:,t)=squeeze(GG(:,:,msel(t)))*ssICt(:,t-1); % Initial conditions - i.e., simulations with no shocks 
    if TT(t)>=2020
        if TT(t)==2020
            ssCt=sst; ssDt=sst;
        end
        % ss6t: fiscal shocks all funded after CARES Act
        ssCt(:,t)=squeeze(GG(:,:,msel(t)))*ssCt(:,t-1)+squeeze(RR(:,:,msel(t)))*SSti(t,size(SSt,2)+1:end)'...
            -1*squeeze(RR(:,3,msel(t)))*SSti(t,size(SSt,2)+3)'+1*squeeze(RR(:,2,msel(t)))*SSti(t,size(SSt,2)+3)'; 
        if TT(t)>=2021 % ss11t: fiscal shocks all funded after ARPA (March 2021)
            ssDt(:,t)=squeeze(GG(:,:,msel(t)))*ssDt(:,t-1)+squeeze(RR(:,:,msel(t)))*SSti(t,size(SSt,2)+1:end)'...
                -1*squeeze(RR(:,3,msel(t)))*SSti(t,size(SSt,2)+3)'+1*squeeze(RR(:,2,msel(t)))*SSti(t,size(SSt,2)+3)';
        elseif TT(t)<2021
            ssDt(:,t)=sst(:,t);
        end
        
    end
end
% Transpose
sst=sst'; ssAt=ssAt'; ssCt=ssCt';ssBt=ssBt'; ssDt=ssDt'; ssICt=ssICt';
if long_sample==1
    %% Save Time Series of Estimated Unfunded Fiscal shocks to be used to simulated the stylized model and obtain Figure X
    filename = '.\Fig_X\EstimatedUnfundedTransferShocks.xlsx';
    xlswrite(filename,[TT',SSti(:,size(SSt,2)+3)]);
    %% Save Time Series of Inflation simuated from the estimated quantitative model using only unfunded fiscal shocks to be used in Figure X
    filename='.\Fig_X\RoleUnfunded.xlsx';
    xlswrite(filename,[TT',4*(ssAt(:,14)-ssICt(:,14))]);
end
if long_sample==1
    %% Figure VII
    
    % Bar Decomposition (inflation)
    R_eb2=[4*(ssAt(:,14)-ssICt(:,14)),4*(sst(:,14))-4*(ssAt(:,14)-ssICt(:,14))-4*(ssBt(:,14)),4*(ssBt(:,14))+model.param(13)*4];
    R_ebneg2=R_eb2;
    R_ebneg2(R_ebneg2>0)=0;
    R_ebpos2=R_eb2;
    R_ebpos2(R_ebpos2<=0)=0;
    
    
    % Bar Decomposition (GDP growth)
    gdp_eb2=4*[(ssAt(:,144)-ssICt(:,144)),(sst(:,144))-(ssAt(:,144)-ssICt(:,144))-ssBt(:,144)-(sst(:,144)-sst(:,69)),(ssBt(:,144))+model.param(12)+(sst(:,144)-sst(:,69))]; 
    gdp_ebneg2=gdp_eb2;
    gdp_ebneg2(gdp_ebneg2>0)=0;
    gdp_ebpos2=gdp_eb2;
    gdp_ebpos2(gdp_ebpos2<=0)=0;
    
    
    temppos_8p = [10, 10, 1050, 1200];
    f=figure('Position', temppos_8p);
    subplot(2,1,1)
    h=bar(TT,R_ebpos2,'stack');
    set(h,{'FaceColor'},{[0,0,0];[.6,.6,.6];[1,1,1]},{'EdgeColor'},{'k';'k';'k'});
    hold on
    h=bar(TT,R_ebneg2,'stack');
    set(h,{'FaceColor'},{[0,0,0];[.6,.6,.6];[1,1,1]},{'EdgeColor'},{'k';'k';'k'});
    hold on
    plot(TT,4*(sst(:,14))+model.param(13)*4,'r','LineWidth',2);
    axis tight
    title('Inflation')
    
   
    subplot(2,1,2)
    h=bar(TT,gdp_ebpos2,'stack');
    set(h,{'FaceColor'},{[0,0,0];[.6,.6,.6];[1,1,1]},{'EdgeColor'},{'k';'k';'k'});
    hold on
    h=bar(TT,gdp_ebneg2,'stack');
    set(h,{'FaceColor'},{[0,0,0];[.6,.6,.6];[1,1,1]},{'EdgeColor'},{'k';'k';'k'});
    hold on
    plot(TT,4*(sst(:,144)+model.param(12)),'r','LineWidth',2)
    axis tight
    h2=legend('Unfunded Transfers Shocks','Other Policy Shocks','Nonpolicy Shocks + Steady State');
    set(h2,'Location','SouthWest')
    title('GDP growth')
    savefigure_pdf(strcat(Figures_path,'\Figure_VII'));
    
    
    %% Figure V
    
    sslin=zeros(length(TT(1:192)),1);
    p1=polyfit([1:61]',sst(1:61,21),1);
    for t=1:61
        sslin(t)=p1(2)+p1(1)*t;
    end
    
    p2=polyfit([62:120]',sst(62:120,21),1);
    for t=62:120
        sslin(t)=p2(2)+p2(1)*t;
    end
    
    p3=polyfit([121:find(TT==2019.75)]',sst(121:find(TT==2019.75),21),1);
    for t=121:find(TT==2022.5)
        sslin(t)=p3(2)+p3(1)*t;
    end
    
    %3-year moving average of the changes in the share of unfunded transfers
    data1=movmean(SSt(2:end,77)+SSt(2:end,22)-SSt(2:end,56)-SSt(1:end-1,77)-SSt(1:end-1,22)+SSt(1:end-1,56),[12,0]);
    %1-year moving average of the ex-ante real interest rate
    data2=4*movmean(sst(1:end,20),[4,0])+log(ssvec1(1))*400-model.param(13)*4;
    
    fig=figure;
    right_color = [1 0 0];
    left_color = [0 0 0];
    
    set(fig,'defaultAxesColorOrder',[left_color; right_color]);
    subplot(1,2,1)
    plot(TT,sst(1:end,21),'-k')
    hold on
    plot(TT,sslin,'--r')
    hold on
    plot([1975.25,1975.25],[-100,150],'--b','LineWidth',1)
    hold on
    plot([1990,1990],[-100,150],'--b','LineWidth',1)
    hold on
    plot([2020,2020],[-100,150],'--b','LineWidth',1)
    hold off
    shadenber
    lgd1=legend('Federal Transfers: Data','Federal Transfers: Trend');
    set(lgd1,'Location','NorthWest','FontSize',10)
    axis tight
    title('Federal Transfers')
    ylabel(['Percentage deviations' newline 'from balanced growth'])
    
    
    subplot(1,2,2);
    yyaxis right
    ha=plot(TT(2:end),data1,'--r');
    hold on
    yyaxis left
    hb=plot(TT,data2,'k');
    hold off
    ylim([-9,15]);
    shadenber
    lgd=legend([ha,hb],{'Changes in Unfunded Transfers (right axis)','Real Interest Rate (left axis)'});

    set(lgd,'Location','NorthEast','FontSize',10)
    
    tmp=corrcoef(data1(1:end),data2(2:end));
    str=sprintf('Correlation= %1.2f',tmp(1,2));
    T = text(min(get(gca, 'xlim')), min(get(gca, 'ylim')), str);
    set(T, 'fontsize', 14, 'verticalalignment', 'bottom', 'horizontalalignment', 'left','FontWeight', 'bold');
    title('Unfunded Transfers and Real Interest Rates')
    ylabel('Percentage points')
    savefigure_pdf(strcat(Figures_path,'\Figure_V'));
    
    
    %% Figure VI
    
    % Share of Transfers: Total transfers | Share of unfunded transfers |
    % Share of funded transfers
    
    z_eb=[SSt(:,76)+SSt(:,22)+SSt(:,77),SSt(:,77)+SSt(:,22)-SSt(:,56),SSt(:,76)+SSt(:,56)];
    
    % Preparing vectors for bar graphs
    z_ebneg=z_eb;
    z_ebneg(z_ebneg>0)=0;
    z_ebpos=z_eb;
    z_ebpos(z_ebpos<=0)=0;
    
    figure
    data3=[SSt(1,77)+SSt(1,22)-SSt(1,56);SSt(2:end,77)+SSt(2:end,22)-SSt(2:end,56)-SSt(1:end-1,77)-SSt(1:end-1,22)+SSt(1:end-1,56)];
    cdata3=cumsum(data3);
    
    % check
    if max(abs([z_eb(:,2)-cumsum(data3)]))>1e-8
        error('gray bars in Figure 5 differ from cumsum(data3)')
    end
    
    subplot(2,4,[1,2,5,6])
    h=bar(TT(1:end),cdata3);
    set(h,{'FaceColor'},{[.4 .4 .4]},{'EdgeColor'},{'k'});
    hold on
    plot([1975.25,1975.25],[0,45],'-.r','LineWidth',1)
    hold on
    plot([1990,1990],[0,45],'-.r','LineWidth',1)
    hold on
    plot([2020,2020],[0,45],'-.r','LineWidth',1)
    hold on
    axis tight
    title('Panel (A): Unfunded Transfers')
    ylabel(['Percentage deviations' newline 'from balanced growth'])
    axes('position',[.345 .78 .12 .12])
    box on % put box around new pair of axes
    indexOfInterest = (TT <= 2022.5) & (TT >= 2019); % range of t near perturbation
    plot(TT(indexOfInterest),cdata3(indexOfInterest),'k'); % plot on new axes
    xticks([2019 2020 2021 2022])
    xticklabels({'2019','2020','2021','2022'})
    axis tight
    grid on
    
    subplot(2,4,3)
    h=bar(TT(1:find(TT==1975)),z_ebpos(1:find(TT==1975),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    h=bar(TT(1:find(TT==1975)),z_ebneg(1:find(TT==1975),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    plot(TT(1:find(TT==1975)),z_eb(1:find(TT==1975),1),'r','LineWidth',1.4)
    axis tight

    
    title('Panel(B) Total Transfers') 
    
    subplot(2,4,4)
    h=bar(TT(find(TT==1975):find(TT==1990)),z_ebpos(find(TT==1975):find(TT==1990),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    h=bar(TT(find(TT==1975):find(TT==1990)),z_ebneg(find(TT==1975):find(TT==1990),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    plot(TT(find(TT==1975):find(TT==1990)),z_eb(find(TT==1975):find(TT==1990),1),'r','LineWidth',1.4)
    axis tight
    title('Panel (C): Total Transfers')
    
    
    
    subplot(2,4,7)
    h=bar(TT(find(TT==1990):find(TT==2019.75)),z_ebpos(find(TT==1990):find(TT==2019.75),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    h=bar(TT(find(TT==1990):find(TT==2019.75)),z_ebneg(find(TT==1990):find(TT==2019.75),2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    plot(TT(find(TT==1990):find(TT==2019.75)),z_eb(find(TT==1990):find(TT==2019.75),1),'r','LineWidth',1.4)
    axis tight
    h1=legend('Unfunded','Funded');
    set(h1,'FontSize',12,'Location','NorthWest')
    title('Panel(D): Total Transfers')
    
    
    subplot(2,4,8)
    h=bar(TT(find(TT==2020):end),z_ebpos(find(TT==2020):end,2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    h=bar(TT(find(TT==2020):end),z_ebneg(find(TT==2020):end,2:3),'stack');
    set(h,{'FaceColor'},{[.4,.4,.4];[1,1,1]},{'EdgeColor'},{'k';'k'});
    hold on
    plot(TT(find(TT==2020):end),z_eb(find(TT==2020):end,1),'r','LineWidth',1.4);
    if long_sample==0
        xticks([2020 2020.25 2020.5 2020.75])
        xticklabels({'2020:Q1','2020:Q2','2020:Q3','2020:Q4'})
    else
        xticks([2020 2021 2022])
        xticklabels({'2020','2021','2022'})
    end
    axis tight
    title('Panel (E): Total Transfers')
    
    savefigure_pdf(strcat(Figures_path,'\Figure_VI'));
    
end
%% Forecasts
nfor=40; ssf=zeros(size(SSt,2),nfor); yyf=zeros(size(ZZ2,1),nfor);ssCf=ssf; yyCf=yyf;
ssDf=ssf; yyDf=yyf;
for i=1:nfor
    ssf(:,i)=GG2^i*SSt(end,:)';
    ssCf(:,i)=GG2^i*ssCt(end,:)';
    ssDf(:,i)=GG2^i*ssDt(end,:)';
end
SSft=[SSt;ssf'];SSCft=[ssCt;ssCf']; SSDft=[ssDt;ssDf'];
for i=1:nobs+nfor
    yyf(:,i)=ZZ2*SSft(i,:)';
    yyCf(:,i)=ZZ2*SSCft(i,:)';
    yyDf(:,i)=ZZ2*SSDft(i,:)';
end


YYft=yyf';
YYCft=yyCf';
YYDft=yyDf';
TTf=[TT(1:end-1),TT(end):.25:TT(end)+nfor/4];

% Transform model's debt into debt to GDP ratio
YY(:,10)=exp((YY(:,10))/100)*25;
YYft(1:length(TT),10)=exp((YYft(1:length(TT),10)+log(model.param(8))*100)/100)*25;
YYCft(1:length(TT),10)=exp((YYCft(1:length(TT),10)+log(model.param(8))*100)/100)*25;
YYDft(1:length(TT),10)=exp((YYDft(1:length(TT),10)+log(model.param(8))*100)/100)*25;


if long_sample==1
    %% Figure VIII
    figure
    subplot(2,2,2)
    plot(TTf(find(TTf==2019):nobs),4*YYft(find(TTf==2019):nobs,7)+2,'k');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*YYCft(find(TTf==2019.75):find(TTf==2022.5),7)+2,'-.r');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*YYDft(find(TTf==2019.75):find(TTf==2022.5),7)+2,'--b');
    hold off
    grid on
    xlim([2019,2022.5]);
    xticks([2019 2020 2021 2022]);
    xticklabels({'2019','2020','2021','2022'});
    title('Inflation')
    
    subplot(2,2,4)
    plot(TTf(find(TTf==2019):nobs),4*SSft(find(TTf==2019):nobs,20)+4*(ssvec1(end)-model.param(13)),'k');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*SSCft(find(TTf==2019.75):find(TTf==2022.5),20)+4*(ssvec1(end)-model.param(13)),'-.r');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*SSDft(find(TTf==2019.75):find(TTf==2022.5),20)+4*(ssvec1(end)-model.param(13)),'--b');
    hold off
    grid on
    xlim([2019,2022.5]);
    xticks([2019 2020 2021 2022]);
    xticklabels({'2019','2020','2021','2022'});
    title('Real Interest Rate')
    
    subplot(2,2,1)
    plot(TTf(find(TTf==2019):nobs),(YYft(find(TTf==2019):nobs,4)),'k');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),YYCft(find(TTf==2019.75):find(TTf==2022.5),4),'-.r');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),YYDft(find(TTf==2019.75):find(TTf==2022.5),4),'--b');
    hold off
    grid on
    xlim([2019,2022.5]);
    xticks([2019 2020 2021 2022]);
    xticklabels({'2019','2020','2021','2022'});
    title('Hours')
    
    subplot(2,2,3)
    plot(TTf(find(TTf==2019):nobs),4*(YYft(find(TTf==2019):nobs,5))+log(ssvec1(1))*400,'k');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*(YYCft(find(TTf==2019.75):find(TTf==2022.5),5))+log(ssvec1(1))*400,'-.r');
    hold on
    plot(TTf(find(TTf==2019.75):nobs),4*(YYDft(find(TTf==2019.75):find(TTf==2022.5),5))+log(ssvec1(1))*400,'--b');
    hold off
    grid on
    xlim([2019,2022.5]);
    xticks([2019 2020 2021 2022]);
    xticklabels({'2019','2020','2021','2022'});
    title('Federal Funds Rate')
    h1=legend('Data','CARES and ARPA fully funded', 'ARPA fully funded');
    set(h1,'FontSize',12,'Location','NorthWest');
    
    savefigure_pdf(strcat(Figures_path,'\Figure_VIII'));
    
    
    
else
    %% Figure IX
    % ARPA Stimulus (2021:Q1)  
    Dataf(end,[2:8,10:22])=NaN; %Take transfer increase (position 9) in 2021 (ARPA stimulus) and NaN all the other observables. 
    QQa=QQ*0; QQa(2,2,1)=QQ(2,2,1);QQa(3,3,1)=QQ(3,3,1); QQa(2,2,2)=QQ(2,2,2);QQa(3,3,2)=QQ(3,3,2); %Shutting down the non-transfers shocks in preparation for the out-of-sample forecasts (only shocks to funded and unfunded transfers allowed to explain the obsetverd growth rate in federal transfers in 2021:Q1 - ARPA shock)
    
    [ SSat,AAat,PPat] = KalmanSmoother( Dataf(end,2:end),GG,RR,QQa,DD,ZZ2,HH,ZZ_t,2,SSt(end,:)',squeeze(PPt(:,:,end))*0 );
    nfor=40;
    for i=1:nfor
        ssaf(:,i)=GG2^i*SSat;
    end
    SSaft=[SSt;SSat';ssaf(:,:)'];
    
    for i=1:nobs+nfor
        yyaf(:,i)=ZZ2*SSaft(i,:)';
    end
    YYaft(:,:)=yyaf(:,:)';
    YYaft(1:length(TT),10)=exp((YYaft(1:length(TT),10)+log(model.param(8))*100)/100)*25;
    
    %%
    %GDP deflator from Fred
    GDPDEFL=[4.9803
        6.2074
        6.0718
        6.6693
        8.1046
        8.8038
        4.2929];
    figure
    subplot(2,2,2)
    
    patch([2020.75 2024.75 2024.75 2020.75],[-2 -2 8.8 8.8 ],'k','EdgeColor','k','FaceAlpha' ,.3,'EdgeAlpha' ,.3,'HandleVisibility','off')
    hold on
    plot(TTf(find(TTf==2019):find(TTf==2020.75)),4*YYft(find(TTf==2019):find(TTf==2020.75),7)+2,'k');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*YYft(find(TTf==2020.75):find(TTf==2024.75),7)+2,'--b');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*YYaft(find(TTf==2020.75):find(TTf==2024.75),7)+2,'-.r');
    hold on
    plot([2020.75:.25:2022.5],[4*YYft(find(TTf==2020.75),7)+2;GDPDEFL],'k')
    

    hold off
    grid on
    axis tight
    h1=legend('Data','Forecast without ARPA','Forecast with ARPA');
    set(h1,'FontSize',8,'Location','SouthEast');
    
    title('Inflation')
    
    subplot(2,2,4)
    patch([2020.75 2024.75 2024.75 2020.75],[-6.9 -6.9 2.8 2.8 ],'k','EdgeColor','k','FaceAlpha' ,.3,'EdgeAlpha' ,.3)
    hold on
    plot(TTf(find(TTf==2019):find(TTf==2020.75)),4*SSft(find(TTf==2019):find(TTf==2020.75),20)+4*(ssvec1(end)-model.param(13)),'k');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*SSft(find(TTf==2020.75):find(TTf==2024.75),20)+4*(ssvec1(end)-model.param(13)),'--b');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*SSaft(find(TTf==2020.75):find(TTf==2024.75),20)+4*(ssvec1(end)-model.param(13)),'-.r');

    hold off
    grid on
    axis tight
    title('Real Interest Rate')
    
    subplot(2,2,1)
    patch([2020.75 2024.75 2024.75 2020.75],[-10 -10 23.5 23.5 ],'k','EdgeColor','k','FaceAlpha' ,.3,'EdgeAlpha' ,.3)
    hold on
    plot(TTf(find(TTf==2019):find(TTf==2020.75)),(YYft(find(TTf==2019):find(TTf==2020.75),4)),'k');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),(YYft(find(TTf==2020.75):find(TTf==2024.75),4)),'--b');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),YYaft(find(TTf==2020.75):find(TTf==2024.75),4),'-.r');

    hold off
    grid on
    axis tight
    title('Hours')
   
    subplot(2,2,3)
    patch([2020.75 2024.75 2024.75 2020.75],[0 0 5 5 ],'k','EdgeColor','k','FaceAlpha' ,.3,'EdgeAlpha' ,.3)
    hold on
    plot(TTf(find(TTf==2019):find(TTf==2020.75)),4*(YYft(find(TTf==2019):find(TTf==2020.75),5))+log(ssvec1(1))*400,'k');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*(YYft(find(TTf==2020.75):find(TTf==2024.75),5))+log(ssvec1(1))*400,'--b');
    hold on
    plot(TTf(find(TTf==2020.75):find(TTf==2024.75)),4*(YYaft(find(TTf==2020.75):find(TTf==2024.75),5))+log(ssvec1(1))*400,'-.r');

    hold off
    grid on
    axis tight
    title('Federal Funds Rate')

    savefigure_pdf(strcat(Figures_path,'\Figure_IX'));
end


